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ABSTRACT 

We first reported evidence for differential rotation of k 1 Ceti in Paper I. In this paper we demon- 
strate that the differential rotation pattern closely matches that for the Sun. This result is based on 
additional MOST (Microvariability & Oscillations of STars) observations in 2004 and 2005, to comple- 
ment the 2003 observations discussed in Paper I. Using StarSpotz, a program developed specifically 
to analyze MOST photometry, we have solved for k, the differential rotation coefficient, and Peq, the 
equatorial rotation period using the light curves from all three years. The spots range in latitude from 
10° to 75° and k — 0.0901q'qq5 - less than the solar value but consistent with the younger age of the 
star, k is also well constrained by the independent spectroscopic estimate of vs'mi. We demonstrate 
independently that the pattern of differential rotation with latitude in fact conforms to solar. 

Details are given of the parallel tempering formalism used in finding the most robust solution which 
gives Peq — 8.77~to;o4 days - smaller than that usually adopted, implying an age < 750 My. Our 
values of Peq and k can explain the range of rotation periods determined by others by spots or 
activity at a variety of latitudes. Historically, Ca II activity seems to occur consistently between 
latitudes 50° and 60° which might indicate a permanent magnetic feature. Knowledge of k and Peq 
are key to understanding the dynamo mechanism and rotation structure in the convective zone as 
well assessing age for solar-type stars. We recently published values of k and Peq for e Eri based 
on MOST photometry and expect to analyse MOST light curves for several more spotted, solar-type 
stars. 

Subject headings: stars: activity - stars: individual: Kappal Ceti, HD 20630 - stars: late-type - 
stars: starspots - stars: rotation - stars: oscillations - stars: exoplanets 



1. INTRODUCTION 

Differential rotation and c onvection provide th e engine 
for the solar dynamo (eg . lOssendrijveil (l2003f0 . The 
mechanism suggested by iLebedinskiil (|1941f) for differ- 
ential rotation seems to be the most widely accepted. 
Angular momentum is continually redistributed in the 
deep convective envelope which occupies some 30% of 
the outer solar radius resulting in nonuniform rotation 
because th e convective turbul ence is affected by the Cori- 
olis force (|Kitchatinovl 120051 ). There is good evidence 
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that a younger Sun would have rotated more rapidly and 
gradually lost angular momentum through magnetic cou- 
pling to the solar wind and there is currently a consid- 
erable theoretical effort to develope quantitative models 
of differential r otation in solar-typ e sta rs (see, for exam- 
ple, re views by iKitchatinovl (|2005f ) and iThompson et al.l 

For stars other than the Sun, the simultaneous detec- 
tion of two or more spots at different latitudes allows a 
direct measurement o f differential rotation . The MOST 
photometric satellite (Walker et al. 2003aj) offers contin- 
uous photometry of target stars with an unprecedented 
precision for weeks at a time making it it possible to 
track s pots accurately . 

We (|Rucinski et al.l 120041 ) have already reported the 
detection of a pair of spots in the 2003 MOST light curve 
of k 1 Ceti (HD 20630, HIP 15457, HR 996; V = 4.83, 
B — V = 0.68, G5V). The photometry was part of the 
satellite commissioning and prior to a considerable im- 
provement pointing in early 2004. The larger spot had a 
well-defined 8.9 d rotation period while the smaller one 
had a period >9.3 d. Only the latitude of the larger spot 
could be deter mined with any confid ence. Spectroscopic 
observations by Shkolni k et alJ (j2003) , mostly from 2002, 
of the Ca II K-line emission rever sals were also synchro - 
nized to a period of about 9.3 d. iRucinski et al.l (|2004f ) 
also determined a stellar radii of i?* = 0.95 ± 0.10 Rq- 

To better decipher the spot activity and rotation of 
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k 1 Ceti and to model spot distributions on other MOST 
targets, one of us, BC, dev e loped t he program Starspotz 
(iCroll et alJ (|2006bD ; iGroi] (|2006aD 1. It was applied first 
to a 2005 MOST light curve covering three rotations of e 
Eri where two spots were detected at different latitudes 
rotating with different periods. From this, we derived a 
differential rotation coefficient for e Eri (see section 4), 
k = 0.II+ Q2, in agreement with models by iBrown et al.l 
( 2004) for young Sun-like stars rotating roughly twice as 
fast as the Sun. 

In this paper, we apply Starspotz to the three separate 
MOST light curves of k 1 Ceti - from 2003, 2004 and 2005. 
We specifically solve for the astrophysically important 
elements of inclination, differential rotation, equatorial 
period and rotation speed and find that the differential 
rotation profile is identical to solar.We use parallel tem- 
pering methods to fully explore the realistic parameter 
space and identify the global x 2 minimum to the data-set 
that signifies a best-fitting unique starspot configuration. 
Markov Chain Monte Carlo (MCMC) methods are then 
used to explore in detail the parameter space around this 
global minimum, to define best-fitting values and their 
uncertainties and to explore possible correlations among 
the fitted and derived parameters. 

2. MOST PHOTOMETRY AND LIGHT CURVES 

The MO ST satellite, launched on June 2003, is fully de- 
scribed bv I Walker et~all (|2003al ). A 15/17.3 cm Rumak- 
Maksutov telescope feeds two CCDs, one for tracking and 
the other for science, through a single, custom, broad- 
band filter (350 - 700 nm). Starlight from primary sci- 
ence targets (V < 6) is projected onto the science CCD 
as a fixed (Fabry) image of the telescope entrance pupil 
covering some 1500 pixels. The experiment was designed 
to detect photometric variations with periods of minutes 
at micro-magnitude precision and does not use compar- 
ison stars or flat-fielding for calibration. There is no di- 
rect connection to any photometric system. Tracking jit- 
ter was dramatically reduced early in 2004 to ~1 arcsec 
which, subsequently, led to significantly higher photo- 
metric precision. 

The observations reported here were reduced by RK. 
Outlying data points generated by poor tracking or cos- 
mic ray hits were removed. At certain orbital phases, 
MOST suffers from parasitic light, mostly Earthshine, 
the amount and phase depending on the stellar coor- 
dinates, spacecraft roll and season of the year. Data 
are recorded for seven Fabry images adjacent to the tar- 
get to track the background. The background signals 
were combined in a mean and subtracted from the tar- 
get photometry. This procedure also corrected for bias, 
dark and background signals and their variations. Re- 
du ctions basically follow ed the scheme already outlined 
bv lRucinski et alJ (|2004f) . 

A log of the MOST observations of k 1 Ceti for each 
of 2003, 2004 and 2005 is given in Table 1. For the 
light curve analysis in this paper, data were binned at 
the MOST orbital period of 101.413 min. Figure 1 dis- 
plays full light curves for each of the three years with 
time given as JD-2451545 (heliocentric). The solid line 



in the lower panel displays the difference between the 
two most divergent background readings from the seven 
background Fabry images. There is no obvious corre- 
lation between structure in the differential background 
signal and the light curve. The complete light curve 
can be downloaded from the MOST Public Archive at 
www.astro.ubc.ca/MOST. 

The formal rms point-to-point precison, a, is given in 
Table 1. The shapes of the light curves ultimately used 
for the StarSpotz analysis depend to a degree on just how 
the background is removed. Despite having the summed 
signal values from 7 adjacent Fabry images we cannot 
interpolate a 'true' background but subtraction of the 
mean background is a consistent option. In each of the 
three data sets used here the 7 background values are 
closely identical when averaged over an orbit but the dif- 
ferences between the extreme backgrounds can be signif- 
icant when the ambient background is high (eg. pointing 
close to the Moon) which is why that difference is shown 
in the lower panel for each light curve. 

We have assumed that the true errors in the individ- 
ual orbital means making up each light curve are pro- 
portional to the difference between the maximum and 
minimum background values. These errors have been 
arbitrarily scaled by one-fourth for the formal StarSpotz 
analysis and these 'errors' are shown as bars in Figures 4, 
5, and 6. The individual bars correspond approximately 
to the standard deviation of the background. 

3. STARSPOTZ 

StarSpotz (jCroll et al.ll2006bl) is a program bas ed on 
SPOTMODEL (Ribari k et al .1120031 : lRibarik|[2002] ). that 
is designed to provide a graphical user interface for pho- 
tometric spot modeling. The functionality included in 
StarSpotz is designed to take advantage of the nearly 
continuous photometry returned by MOST to quantify 
the non-uniqueness problem often associated with pho- 
tometric spot mod eling. B a sically , S tarSpotz uses the 
analytic models of iBuddinel (| 19771 ) or iDorrenl (|1987l ) to 
model the drop in intensity caused by Nspqt circular, 
non-ov erlapping spots. The model proposed in iBuddina 
(|1977|) is used in this application. 

Fitting can be performed by a variety of methods in- 
cluding a Marquardt-Levenberg non-linear least squares 
algorithm (Marquardt 1963; Levenberg 1944; Press et al. 
[1992) , or Markov Chain Monte Carlo (MCMC) function- 
ajitv (|Gregorvll2005bl ; iGamerman et alj|1999t iGilks et all 
1996). The parameters of interest are the stellar incli- 
nation angle, i, the unspotted intensity of the star, U , 
and the period, epoch, latitud e and angular size of the 
a th spot: p a , E a , /3 a , and j a . ICroll et all (|2006b| ) have 
demonstrated how the Marquardt-Levenberg non-linear 
least squares algorithm c ould be used to fit the cycle to 
cycle variation in e Eri. Croll ( 2006a]) investigated the 
same observations of e Eri demonstrating that Markov 
Chain Monte Carlo functionality is very helpful in defin- 
ing correlations and possible degeneracies between the 
fitted parameters. 

The search for additional local minima in y 2 space indi- 
cating other plausible spot configurations that produce a 
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TABLE 1 
MOST Observations of k 1 Ceti 



year 


dates 


total 


duty cycle 


exp 


cycle 




start 


HJD - 2451545 


days 


% 


sec 


sec 


mmag 


2003 5 Nov 


1403.461 to 1433.901 


30.44 


96 


40 


50 


0.5 


2004 15 Oct 


1748.578 to 1768.989 


20.41 


99 


30 


35 


0.30 


2005 14 Oct 


2112.721 to 2125.065 


12.34 


83 


30 


35 


0.19 



for data binned at the MOST orbital period of 101.413 min. 



reasonable fit to the light-intensity curve is an important 
part of the modeling because it allows one to quantify 
the uniqueness of a solution. We argue that this task 
of completely searching parameter space for other phys- 
ically plausible local minima i s better suited t o a pro- 
cess ca lled parallel tempering ()Gregorvl l2005al ; I Gregory! 
I2005bf) . discussed he re in jj4.2l rather th an the Unique- 



ness test methods in lCroll et al.l ( 2006b). 

The present and subsequent releases of the freely avail- 
able program StarSpotz, including the full source-code, 
executable, and documentation are available at: 
www.astro.ubc.ca/MOST/StarSpotz.html 
Release 3.0 includes the functionality discussed which 
allows one to fit multiple epochs of data as well as the 
parallel tempering functionality. 

4. THE SPOTS ON k 1 CETI 

As the MOST observations of k 1 Ceti come from three 
well-seperated epochs with independent starspot groups 
from epoch to epoch, the fitting process w as consider- 
ably mo re c omplicated th an those detailed in lCroll et al.l 
(|2006bJ ) and lCrollI (j2006af ) . Rather than fitting for the pe- 
riod, p, of the individual spots in each epoch we derive a 
common equatorial period, Peq, and the differential ro- 
tation coefficient, k, among the three epochs where these 
are given by: 

Pp = Peq/(! " k sin 2 /?) (1) 
where Pp is the period at latitude f3. There are five fewer 
parameters to be fitted because the periods of the in- 
dividual spots are no longer independent variables, but 
completely defined by k, Peq, and f3. Analyses of so- 
lar differential rot ation usually include terms in sin 4 /3 
( Kitchatinov 2005) but we did not expect to detect spots 
at such high latitudes where this term would be impor- 
tant. 

Common limb-darkening coefficients, u = 0.684, and 
flux ratios between the spot and unspotted photosphere, 
Ku, = 0.22 were ass umed. The former was b ased on 
the compilations of iDiaz-Cordoves et alj (|1995l ). while 
the latter comes from the canonical value for sunspots 
(T spof ~ 4000 K, T phot ~ 5800 K). The 27 parameters to 
be fitted are summarized in Table [2l 

Initial investigations of the light curves using the 
StarSpotz standard spot model indicated that different 
numbers of spots were required for the various epochs. 
The 2003 and 2005 MOST data-sets could be well ex- 
plained by two starspots rotating with different periods, 
at different latitudes. The behaviour observed in the 
2004 MOST light curve, on the other hand, required a 



third spot. Thus the following analysis explicitly assumes 
two, three and two starspots for the 2003-2005 MOST 
epochs, respectively. 

In this present application we use the same para digm 
that wa s applied to re 1 Ceti bvlRucinski et al.1 (l2004f ). and 
e Eri bv lCrol (pOOpah and lCroll et all (|2006bh where the 
spot parameters (E, (3, 7) were assumed to be constant 
throughout the observations, and that the cycle-to-cycle 
variability observed is due solely to spots rotating dif- 
ferentially. Other interpretations are certainly plausible 
including circular spots with time- varying spot sizes and 
latitud es such as those applied to IM Peg in Ribarik et al. 
(2003). The Maxi mum Entropy fc Tik honov reconstruc- 
tion technique of lLanza et al.l (|1998[ L which assumes 
pixellated spots of various contrasts, is another example. 
We feel that with the limitations of photometric spot 
modeling, our assumption of unchanging spot parame- 
ters is the most suitable for the timescale of the MOST 
light curves. It provides an opportunity to determine 
a unique solution, as well as recovering astrophysically 
important parameters. 

We recognise, given the obvious changes in the k 1 Ceti 
light curve over a year, that the assumption of spot 
parameters constancy over a few weeks is a simplifica- 
tion but we do not expect significant spot migration on 
timescales of less than a month. 

4.1. StarSpotz Markov Chain Monte Carlo (MCMC) 
functionality 

The Markov Chain Monte Carlo (MC MC) fu n ctiona l- 
ity used in StarSpotz is fully described in ICrollI (|2006af h 
The goal is to take an n-step intelligent random walk 
around the parameter space of interest while recording 
the point in parameter space for each of the K fitted 
parameters for each step. This samples the posterior 
parameter distribution, which can be thought of in pho- 
tometric spot modeling as simply the range in each of the 
K fitted parameters that provide a reasonable fit to the 
light curve. The advantage of a MCMC is that it does not 
simply follow the path of steepest descent - resulting in 
the possibility that it becomes stuck in a local minimum 
- but allows the MCMC to fully explore the immediate 
parameter space. The K parameters can thus be defined 
by a K dimensional vector y. 

T stands for the "virtual temperature" and determines 
the probability that the MCMC chain will accept large 
devia tions to higher regions in \ 2 space ()Saiina et al.l 
2006) . The virtual temperature should properly be unity 
(T = 1) to sample the posterior parameter distribution 
when reduced x 2 is near unity. However, in photomet- 
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ric spot modeling, nature is often more complicated than 
one's model can take into account, resulting in anoma- 
lously high reduced x 2 values. Setting T to be the mini- 
mum reduced x 2 value observed produces the same effect 
as scaling reduced \ 2 to one, and thus effectively sam- 
ples the posterior parameter distribution. In this present 
application T m ; n will refer to this value of T equal to the 
minimum reduced \ 2 observed. 
A suitable burn- in period as described in lCrollI ( 



is often used, or a x cu t where points with \ 
values above this user-specified minimum are excised 
from the ana lysis. A user-specified thinning factor, /, 
(|Crolll I2006ah is also used. To ensure proper conver- 
gence and mixing we em ploy the tests as proposed by 
iGelman and Rubinl (Il992f) and explicitly defined for our 
purposes in ICrofF(|2006a| ). We ensure the vector R 
falls as close to 1.0 as possible for each of the K fitted 
parameters. The marginali zed likelihood (ISaiina et al. 
l2006HLewis and Bridlell2002l ) as described in §3.1 oflCrjjl 
(|2006a| ) is used to de fine our 68% credi ble regions, while 
the mean likelihood (|Saiina et al.ll2006f ) is used for com- 
parison. The mean likelihood is produced using the % 2 
values in a particular bin of a histogram of parameter 
values, while the marginalized likelihood directly reflects 
the fraction of times that the MCMC chains are in a 
particular bin. Peaks in the marginalized likelihood that 
are inconsistent with peaks in the mean likelihood can be 
indicative of the chains having yet to converge, or larger 
mu ltidimensional paramet er space rather than a better 
fit (|Lewis and B ridle 20Q2) ■ In our analysis we quote the 
minimum to maximum range of the 68% credible regions. 

4.2. Parallel Tempering 

iGregorvi (|2005af ) applied parallel tempering to a radial 
velocity search for extrasolar planets, as the parameter 
space he investigated had widely seperated regions in pa- 
rameter space that provided a reasonable fit to the data. 
In our application as we fit for a large number of param- 
eters, K=27, we use parallel tempering to explore our 
multi-dimensional parameter space and search for other 
local minima in x 2 space indicative of other plausible so- 
lutions, and thus starspot configurations, that produce a 
reasonable fit to the light-intensity curve. 

In parallel tempering, Mp multiple MCMC chains are 
run simultaneously each with different values of the vir- 
tual temperature, T m , where m is an index that runs 
from 1 to Mp. Chains with higher values of T are 
easily able to jump to radically different areas in pa- 
rameter space ensuring that the chain does not become 
stuck in a local minimum. Chains with lower values of 
T can refine themselves and descend into local minima 
and p ossibly the global minimum in x 2 space ((Gregory! 
I2005al) . The chain with T = 1.0xT min is the target dis- 
tribution and is referred to as the "cold sampler." The 
parameters of the "cold sampler" are used for analysis 
to return the posterior parameter distribution. A set 
mean fraction of the time, Z, a proposal is made for 
the parameters y of adjacent chains to be exchanged. 
If this proposal is accepted than adjacent chains are cho- 
sen at random and the K parameters y for each chain 



are exchanged. In this way, this method allows radi- 
cally different areas in parameter space to be investi- 
gated in the high temperature chains, before reaching 
the "cold sampler" where a local minimum can be sought 
efficiently. We run 11 parallel chains with virtual temper- 
atures, T m , of: T = {100xT mm , 50.0xT mm , 33.0xT„ lm 
, 20.0 x Twin, lO-OxT^^, 5.0 %T m i n , 3.00 xT m j n , 
2.00xX , m ^ n , 1.50x J m £ n 1.25 x J m i n . 1.00 xT min }. 

In a typical parallel tempering scheme the number of 
iterations between swap proposals is set to some con- 
stant integer, B (e.g. B w 10). With the large number 
of parameters being fitted across the three epochs of data 
(K=27) in the current application, it was found that this 
typical scheme was not able to efficiently reach the low-% 2 
valleys in the multi-dimensional parameter space indica- 
tive of a reasonable fit to the light-intensity curve. Thus 
the number of iterations between swap proposals, B m , 
was set to a variable number between the chains. B rn 
was set to a larger number for the lower temperature 
chains to allow these chains a greater number of itera- 
tions to refine themselves and explore the low-% 2 areas. 
This was found to be a suitable compromise as the higher 
virtual temperature chains required fewer iterations to 
reach radically different areas of parameter space, while 
the lower virtual temperature chains required nearly an 
order of magnitude more iterations to reach lower ar- 
eas in x 2 space indicative of reasonable fits to the light- 
curve. The number of chain points per chain was set 
to B {B m i n: B m i n: B m i n , B Ta j jTll B Ta j jTll B m i n: B m i n: 
2.0xi? m j n , 2.0xB m i ni 3.0xB m i n , 4. Ox B m i n } . B m i n was 
set to twenty (B m j n =20 ), within an ord er of magnitude 
of the value proposed by I Gregory! (|2005afk extensive tests 
have proven this choice to be reasonable. / was set as 0.80 
to allow for exchanges between adjacent MCMC chains 
on average 80% percent of the time. That means a swap 
is approved if ui < I, where ui is a random number 
(ui e[0, 1]). Significant experimentation was needed to 
choose and refine suitable choices of T m , B m , B m i n , and 
I. 

In the present application we start 8x11 of these par- 
allel tempering chains starting from random points in 
our acceptable parameter space as given in the Paral- 
lel Tempering initial allowed range column of Table O 
MCMC fitting was then implemented for n=4000 steps 
to ensure that the starting point for the parallel tem- 
pering chains provided a mediocre fit to the light curve. 
A mediocre fit was desired to ensure that the starting 
point for the parallel tempering chains provided a rea- 
sonable fit to the light curve so as to reduce the required 
burn-in period, while not starting from within a signifi- 
cant local minimum so as to allow the parallel tempering 
chains to efficiently explore the parameter space. Each of 
these 8x11 MCMC parallel tempering chains were then 
run for n.^13000 exchanges, resulting in n K x20x4 = 
1.04xl0 6 total steps for each of the 8 "cold samplers." 
A burn-in period was not used in this application for our 
analysis. Rather a x 2 cu t was used where all points of 
the parallel tempering "cold samplers" with reduced x 2 
greater than some particular minimum, Xcut> were ex- 
cised from the analysis. This serves a similar function to 
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a burn-in period for the initial points, while allowing the 
analysis to focus on the chain points that provide a rea- 
sonable fit to the light curve. Given the large number of 
parameters we fit for, if =27, this is useful as the paral- 
lel tempering chains spend an inordinate amount of time 
exploring parameter spaces that do not provide a good 
fit to the data, and would otherwise complicate the anal- 
ysis. A thinning factor, /, of 20 were used in all chains 
for the parallel tempering analysis. The 8x11 parallel 
tempering chains required a total computational time of 
7.0 CPU days on a Pentium processor with a clock speed 
of 3.2 Ghz with 1.0 GB of memory. This value of n x was 
chosen to ensure that the vector R, in t his case Rparaiiei 
(jGelman and Rubinl [199^ : ICrollll2006al) . fell as close to 
1.0 as possible for each of the K = 27 parameters. R 
< wl.3 signifies that the chain points have roughly con- 
verged and thus we are close to an equilibrium distribu- 
tion. To determine Rparaiiei we do not use the above x 2 
cut but instead use a burn-in period of ni, urn = 30000 
steps. This is because we are interested in ensuring the 
parallel tempering chains have converged while exploring 
our parameter space as a whole, rather than the narrow 
region of parameter space signifying the global x 2 mini- 
mum. Analysis of these 8 "cold samplers" is given below 

in gjro 

4.2.1. k1 Ceti Parallel Tempering results 

The Parallel Tempering results as applied to MOSTs 
2003-2005 observations of k 1 Ceti indicated that there is 
a single unique solution to the data-set that provides the 
best fit. Following the 71^=13000 exchanges as described 
above the vector Rp ara iiei was investigated for each of 
the K=27 parameters. It was found that Rparaiiei fell 
below 1.6 for all parameters, and Rparaiiei fell below 1.3 
for most parameters. These low values of Rparaiiei indi- 
cate adequate convergence has been obtained. Each of 
the components of the vector Rparaiiei are given in Table 
[2j These low values for all the K parameters of Rparaiiei 
indicate that our parallel tempering chains have approxi- 
mately converged to an equilibrium, and adequately sam- 
pled the realistic parameter space. 

For our analysis we set T m i n as 2.96 as it is the global 
minimum observed in the below section ( §4.3| and a re- 
duced Xcut °f 3-38 as it was judged this was reasonably 
close to the global x 2 minimum to focus the analysis ex- 
clusively on the solutions that provided a reasonable fit 
to the light curve. Our 68% credible regions are given 
by the marginalized likelihood following the x 2 cu t. The 
immediate parameter space indicated by these results are 
given in Table [3] This parameter space was thus worthy 
of more detailed investigation to explore this x 2 global 
minimum, as well as to define best-fit values and appro- 
priate uncertainties for the fitted and derived parameters 
of interest. Also, possible correlations between the fitted 
parameters were investigated for. 

Detailed investigation of the parallel tempering results 
indicate that all other local minima produce significantly 
worse fits to the data-sets. The only minor exception 
indicated by the parallel tempering results is that the 
second spot in the 2003 epoch can also be found at a 



latitude, /?2003_2, hi the southern hemisphere as well as 
the northern hemisphere. Detailed investigation indicate 
that the northern hemisphere solution provides a signif- 
icantly better fit and thus is the solution that will be 
investigated. Given the fact that all other local minima 
have been ruled out, and the low values of Rparaiiei for 
all K parameters, we feel fully justified in stating that 
this parameter space is the unique realistic solution to 
the observed light curve given our assumptions. 

4.3. k 1 Ceti MCMC application 

In iJ4.2l we demonstrated that the solution presented 
is unique for the 2003-2005 MOST epochs under the k, 
Peq, par adigm Thus we used the MCMC methods dis- 
cussed in lCrollI (|2006af ) and summarized above in ij4.1l to 
investigate possible correlations and derive best-fit values 
and uncertainties in the fitted and derived parameters for 
k 1 Ceti. We used M=4 chains starting from reasonably 
different points in parameter space, although within the 
local x 2 minimum that defined the parameter space of 
interest of H4.2I The prior ranges are given in Table [H 
The choices in prior are identical to those used for the 
Parallel Tempering section except a slightly more con- 
servative range is used for the inclination angle, 30° < i 
< 80°, to limit the MCMC chain to a realistic parameter 
space. We run each of the M=4 chains for n=7.0xl0 6 
steps, requiring a total computational time of 8.0 CPU 
days on the aforementioned processor. This large value of 
n was motivated to ensure that the R vector, in this case 
Rmcmc, fell below 1.20 for all K=27 parameters, and 
fell below 1.05 for most parameters, indicative of suitable 
convergence. The values of each of the K components of 
Rmcmc are given in Table [H 

The results of this analysis, henceforth referred to as 
Solution 1, can be seen in Figures 2 and 3, and are sum- 
marized in Tabled As one m ight expect, e specially given 
the similar results of e Eri (jCrollI [20 06a). the latitudes 
of the various spots are moderately correlated with the 
inclination angle of the system, i. However, the differen- 
tial rotation coefficient, fc, is largely uncorrelated with i. 
This is in stark contrast to the results of e Eri where the 
correlation of spot latitude with i produced a moderate 
correlation between k and i. Our derived value of k = 
0.085 - 0.096 is therefore very robust as it is does not 
require an independent estimate of i. 

It is important to note, though, that the 2003 and 2005 
MOST data-sets do not severely limit the allowed best- 
fit range of the differential rotation coefficient, k. The 
2003 and 2005 MOST data-sets were best-fit with spots 
that were not widely seperated in latitude, and thus these 
epochs contributed only marginally to the determination 
of the differential rotation coefficient. The 2004 MOST 
data-set most severely limited the best-fit range of the 
differential rotation coefficient, as it required three spots 
ranging in latitude from the mid-southern hemisphere, 
to the equator, to near the north pole. 

Also, although we have used a relatively liberal flat 
prior on the stellar inclination angle, i, of 30° < i < 
80°, the value of i returned from our fit gives wsini = 
4.64 - 4.94 km s , entirely consistent with the spectro- 
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TABLE 2 
k 1 Ceti MCMC fitted parameters 



parameter 


prior allowed 
range 


Parallel Tempering 
initial allowed range 


Rparallel 


Rmcmc 


1 ( ) 


nr orz on on a, 

25 - 80, 30 - 80 


25 - 85 


1.033 


1.010 


k 


n "7r n 

— 0.7o - U.75 


—0.6 - 0.6 


1.443 


1.053 


F EQ (a) 


en i n k 
o.U - 1U.0 


o.U - 1U.0 


1.4/ / 


1.U1Z 


T J 


n on i m 


1.000 


1 A 1 K 

1.415 


1.024 




1 /i nn n i /ii q n 
14Uy.U - 14io.U 


1 a nn n 1 A i q n 
14Uy.U - 141o.U 


l.OOO 


l.UUD 


P2003_l ( ) 


—90 - 90 


— 10. U - bO.O 


1.144 


1.023 


72003-1 { ) 


U.U - oU 


4.U - o.U 


i 1 1 n 
1.11U 


l.Uol 


-^2003.2 (JD-2451545) 


1 ac\a f\ i /ino n 
14U4.U - 14Uo.U 


i ac\a n i/inon 
14U4.U - 14Uo.U 


1.051 


1.009 


P2003.2 ( ] 


— 90 - 90 


inn £>n n 

— 1U.U - uU.U 


1.548 


1.026 


72003.2 (°) 


U.U - oU 


/in en 
4.U - o.U 


1 Qli) 

l.olz 


i no*; 
l.Uzo 


t/2004 


n on i m rr 

u.yy - 1.015 


1.000 


1.318 


1.124 


#2004.1 (JD-2451545) 


1754.7 - 1758.7 


1754.7 - 1758.7 


1.179 


1.004 


P2004.1 ( ) 


—90 - 90 


1 n n an n 
— 1U.0 - 6U.0 


1.338 


1.013 


72004.1 (°) 


0.0 - 30 


2.5 - 6.0 


1.126 


1.022 


#2004.2 (JD-2451545) 


1761.5 - 1765.5 


1761.5 - 1765.5 


1.391 


1.005 


P2004.2 (°) 


-90 - 90 


-10.0 - 60.0 


1.596 


1.098 


72004.2 (°) 


0.0 - 30 


2.5 - 6.0 


1.209 


1.100 


#2004.3 (JD-2451545) 


1767.1 - 1771.1 


1767.1 - 1771.1 


1.079 


1.039 


/92004.3 (°) 


-90 - 90 


-10.0 - 60.0 


1.490 


1.042 


72004.3 (°) 


0.0 - 30 


2.5 - 6.0 


1.237 


1.035 


t%005 


0.99 - 1.01 


1.000 


1.235 


1.179 


£2005.1 (JD-2451545) 


2115.0 - 2119.0 


2115.0 - 2119.0 


1.182 


1.057 


/32O05.1 (°) 


-90 - 90 


-10.0 - 60.0 


1.249 


1.059 


72005.1 (°) 


0.0 - 30 


4.0 - 8.0 


1.122 


1.070 


#2005.2 (JD-2451545) 


2118.1 - 2122.1 


2118.1 - 2122.1 


1.289 


1.038 


P2005.2 (°) 


-90 - 90 


-10.0 - 60.0 


1.074 


1.072 


72005.2 (°) 


0.0 - 30 


4.0 - 8.0 


1.057 


1.092 


a The left hand values 
application f £|4.3|I. 


were used for the parallel tempering 


application (34.2.10. while the values 


on the right were 


used for the MCMC 



scopic value oiysmi = 4.64 ± 0.11 km s determined by 
IRucinski et~al] (|2004h from the width of the Ca II K-line 
reversals. This agreement is a very useful sanity check 
and gives us great confidence that our fitted results and 
our assumptions of circular spots are in close agreement 
with the actual behaviour of the star during these three 
epochs. 

The minimum \ 2 point observed in our MCMC chains 
is also quoted in Table [3] The fit to the light curves in 
2003, 2004, 2005 and the views of the modelled spots 
shown from the line of sight (LOS) of this minimum \ 2 
solution are given in Figures 111 [5l and [5] 

Our current results are largel y consistent with the an al- 
ysis of the 2003 MOST data bv IRucinski et all ([20041 ) as 
can be seen in the comparison given in Table [3] The 
only discrepancies of note are that our MCMC analysis 
indicates a slightly smaller value of the inclination an- 
gle, i = 57.8 — 63.5°, and that we find a shorter period 
of the 2 n d spot, P2003-2 = 9-02 2 - 9.094 d. The period 
found bv IRucinski et alj (|2004l) for the 2 nd spot depends 
on effective removal of the variations caused by the larger 
spot and it was not possible to asssign a reliable latitude. 
They also asumed that the spots were black making them 
smaller than those in this paper. 
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4.4. k 1 Ceti 2004 Parallel Tempering and MCMC 
application 

Given the good agreement of the MOST 2003, 2004 
and 2005 n 1 Ceti observations with the assumed solar- 
type differential rotation profile in Equation 1 and sum- 
marized in Figure [71 we decided to test independently 
whether solar-type differential rotation is indeed present 
in k 1 Ceti. The 2004 light curve was chosen because, as 
noted above, it constrained the equatorial period, Peq, 
and differential rotation coefficient, k, more significantly 
than the light curves from the other two years. 

For this independent test we fitted explicitly for the 
periods of each of the three spots: P2004_ij P2004_2, and 
P2004_3- The K=13 fitted parameters are summarized in 
Table [H For these we used the same priors, and for the 
parallel tempering chains we started from the same ran- 
dom ranges in parameter space, as summarized in Table 
O We assumed i= 60.6° corresponding to the minimum 
X 2 value found above. This assumption for the inclina- 
tion angle is justified because it results in wsini = 4.77km 
s _1 (assuming P E q rj 8.7 8), a value close to the spectro- 
scopic value measured bv lRucinski et al.l (|2004[ ). 

Parallel tempering as described in N4.2I was used to 
find the unique global minimum. Starting from random 
points in acceptable parameter space we implemented 
MCMC fitting for n=3000 steps to ensure a mediocre 
fit to the light curve. Parallel tempering chains were 
then implemented with n x = 1800 exchanges resulting in 
144000 steps for each of the 8 "cold samplers" . We used 
a burn-in period of the first 5000 steps to determine the 
R vector, in this case Parallel- All K=13 parameters 
of Rp'araiiei f en below 1.4, and below 1.2 for most, as 
summarized in Table 0J Analysis is performed using a 
T m in of 1.50 and a reduced Xcut 01 2.0. The parallel 
tempering results identify a unique global minimum as 
summarized in Table [H that is nearly identical to Solu- 
tion 1 as summarized above in t j4.3l and Table 3. 

Thus after using parallel tempering to identify the 
unique global minimum, we explored this parameter 
space with our MCMC techniques as described above 
in 34. II and M4.3I We use n=3.7xl0 6 steps, and then use 
a burn-in period of 1000 steps to determine the R vec- 
tor (in this case P^mcmc)- The parameters of R'mcmc 
fell below 1.02 for all K—13 parameters as summarized 
in Table 01 and thus demonstrated more than adequate 
convergence. Our MCMC results are also summarized in 
Table E] 

As summarized in the bottom-panel of Figure [7] these 
results indicate that the 2004 MOST data-set provides a 
strong independent argument that solar-type differential 
rotation pattern defined by Equation 1 has been observed 
in k 1 Ceti. Using only the 2004 data-set the values of 
the differential rotation coefficient, and equatorial period 
are: k = 0.092 - 0.101, and Peq = 8.65 - 8.72d. These 
values were determined by averaging for each point of 
the MCMC chains the three possible k and Peq values, 
between spots 1 & 2, 1 & 3, and 2 & 3. 

5. DISCUSSION 

The MCMC Solution 1 in Table 3 provides our best 
values for the various spot and stellar parameters. The 



solutions are expressed by the 68% credible ranges rather 
than best values with standard deviations since the like- 
lihood histograms are non-gaussian. Figure [7J is a good 
summary of our analysis. The periods and 68% 
marginalized contours are shown for each of the spots to- 
gether with mean and limiting curves for k (0.090, 0.085, 
0.096) and the range of Peq (8.77, 8.74, 8.81d). Each 
of the three years is distinguished by a different color. 
|/3 1 ranges from 10° and 75° with the 2004 data obvi- 
ously providing the most r igorous constra i nt on k. In- 
deed, it demonstrates that [Rucinski litTaTj (|2004l ) had a 
most challenging task in demonstrating differential ro- 
tation because the spots in 2003 were at very similar 
latitudes with the second spot being particularly small. 

The analysis defined k using Equation 1 which is de- 
rived from the solar pattern. The agreement of all seven 
spots with the form of the k curve as well as the good 
agreement of the spot solutions in Figures 5, 6 and 7 sug- 
gests that the the differential rotation curve for k 1 Cet 
is closely similar to solar. The analysis of the 2004 light 
curve independently of any assumption about k offers 
strong confirmation that the pattern is indeed solar. 

For the Sun, k has been derived quite inde- 
pendently from either su nspot latitudes and periods 
(iNewton and Nunnl 1 1 9 5 lh or surface radial velocities 
([Howard et al.l[l983l) yield ing k = 0.19 and 0.12, respc- 
tively ([Kitchatinovl I2005T ). The large discrepancy be- 
tween these values may be related in part to the different 
behavior of sunspot and photospheric motions. In our 
case, we depend on spots to determine k and our value 
of 0.09 is significantly lower than either of those for the 
Sun. This is in line with calculations by iBrown et al.l 
(2004 ^ who find that k should increase with age. 

iGudel et all (|1997h estimated an age of 750 Myr for 
k 1 Cet from their estimated 9.2 d rotation period. Our 
value of Peq — 8.77 d suggests a still younger age. 

All of the photometric periods found to date for k 1 Cet 
can be explained by spots appearing at different lati- 
tudes . A period of 9.0 9 d is given in the Hipparc os cat- 
alog (|ESAlll997h w hile iMessina & Guinanl (|2002f > quote 
a value of 9.214 d. Bal iunas et a l. (1995) monitored Ca 
II H & K photoelectrically between 196 7 and 1991 and 
found a ro tational period o f 9.4 ± 0.1 d (Bal iunas et al.l 
Il983h and IShkolnik etafl (j2003l ) found a closely simi- 
lar period of ~9.3 d mostly from spectra taken in 2002. 
While the apparent persistance of this period which cor- 
responds to a range of 50° to 60° in latitude for some 35 
years maybe fortuitous, it might be related to some large 
scale magnetic structure. 
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TABLE 3 
re 1 Ceti MCMC fitted parameters 



parameter 


111 ICU 


„ ... 

JrvUClIlSKl 


TJ 11 1 T 

r^arallel J. empermg 


... . 2 

IVlimmum v 








ZUU^i I eo Lll La 


results 


solution 


Solution 1 


Reduced x 2 


n/a 




3.145 - 3.374 


2.936 


2.968- 3.018 


v a 


n/a 


_ 


817 


817 


817 


i (°) 


yes 


70 ±4 


44.9 - 74.1 


60.6 


57.8 - 63.5 


k 


yes 




0.094 - 0.117 


0.087 


0.085 - 0.096 


Peq 


yes 


_ 


8.71 - 8.90 


8.784 


8.74 - 8.81 


vsixii (km s — ^) c 


derived 


4.64 ± 0.11 b 


3.83 - 5.34 


4.77 


4.64 - 4.94 


equatorial speed (km s — ) c 


derived 




5.40 - 5.53 


5.47 


5.46 - 5.50 


u 


assumed 


0.80 


0.6840 


0.684 


0.684 


Km 


assumed 


0.00 


0.220 


0.220 


0.220 




VPS 




0.9964 - 1.0031 


1.0011 


1.0003 - 1.0017 


&?nm 1 (JD-2451545) 






1410.94 - 1410.97 


1410.93 


1410.92 - 1410.94 


£>2003_i (days) 


derived 


s q 4- n 1 d 
o.y in u. 1 


Q Q7Q Q 007 


y.uuo 


q nqn q m ^ 
o.yyu - y.uiij 


P2003_l I ) 


yes 


An 4- 7 


1 ft 7 Qft Q 

lu. / - ou.y 


32 4 


on t; oa 


72003-1 (°) 


yes 


11.1 i 0.6 


1 1 n^i 1 1 qs 


11.75 


1 1 fiS 1 1 Kfi 

1J-.UO ll.OU 


J - / 2()()3_2 yo i-J- £*±a ±0*±0 J 


yes 




1 Af)^ Rt\ 1 zL0.fi 9S 


1 40fi 1 n 




P2003_2 (days) 


derived 


n q o t d 
y.o — y. ( 


Q O/l Q 01 ftn 


y.u / U 


y.uzz - y.uy4 


to \ 

P2003_2 ( ) 


yes 




-45.7 - 00. 9 


37.2 


00 n or\ 

32.9 - 39.8 


72003-2 ( ) 


yes 




no ,4 1 
0.U0 - y.41 


o.yo 


^ ftQ ft 1 Q 
O.OO - O.lO 


T J 

U200A 


yes 


n/a 


1 nnoc; 1 nnno 

i.uuzo - l.uuyz 


1.0149 


1 m on 1 m c:n 
i.uizy - l.UloU 


-^2004_1 (,J-L»-^40l040j 


yes 


n/a 


1 v^ft fto 1 v^ft 7C 

l/OD.Dy - IfOO. (0 


1 / 00. Do 


1 v^ft ftft 1 v^ft 7n 


P2004_i (days) 


derived 


n/a 


701 non 

0. (21 - c.yoy 


8.820 


o.7o7 - o.ool 


P2004_l { ) 


yes 


n/a 


n e 1 c ri 
-U.o - 10. U 




n q 1 ft q 

y.o - lo.o 


72004.1 ( ) 


yes 


n/a 


/ .1 ( - O.Ol 


1 no 


1 10 q no 

/ . ( - o.uy 


"2004.2 (, J JJ- Z401 040 J 


yes 


n/a 


1 "7ft Q c;Q 17fiO ftE 
1 (OO.Oo - 1 (DO. OO 


1 ( 00.04 


1 /in 1 'zftQ c;ft 
1/Dv).4y - 1 /OO.OO 


P2004_2 ^wtyo^ 






q 0S9 - Q SI fi 


9.200 


Q IflQ Q 9S1 

C7.J_iJO C».ivi_>J_ 


#2004.2 (°) 


yes 


n/a 


-43.9 - -28.7 


-46.3 


-47.6 - -43.2 


72004.2 (°) 


yes 


n/a 


7.14 - 19.21 


16.76 


14.44 - 17.31 


#2004.3 (JD-2451545) 


yes 


n/a 


1769.01 - 1769.32 


1769.09 


1769.03 - 1769.16 


P2004.3 (days) 


derived 


n/a 


9.572 - 9.751 


9.572 


9.542 - 9.640 


#2004.3 (°) 


yes 


n/a 


59.2 - 76.9 


77.3 


74.9 - 78.4 


72004.3 (°) 


yes 


n/a 


7.53 - 13.23 


13.04 


11.60 - 13.53 


^2005 


yes 


n/a 


1.0000 - 1.0053 


1.0042 


1.0029 - 1.0051 


E2005.1 (JD-2451545) 


yes 


n/a 


2116.96 - 2117.09 


2116.94 


2116.95 - 2117.05 


P2005.1 (days) 


derived 


n/a 


9.340 - 9.419 


9.372 


9.365 - 9.423 


#2005.1 (°) 


yes 


n/a 


41.8 - 63.0 


58.4 


55.4 - 62.1 


72005.1 (°) 


yes 


n/a 


8.30 - 9.37 


9.74 


9.24 - 10.28 


&2005.2 (JD-2451545) 


yes 


n/a 


2120.09 - 2120.21 


2120.22 


2120.21 - 2120.30 


P2005.2 (days) 


derived 


n/a 


9.094 - 9.253 


9.248 


9.166 - 9.252 


#2005.2 (") 


yes 


n/a 


28.0 - 56.0 


49.6 


42.9 - 50.1 


72005.2 (°) 


yes 


n/a 


7.69 - 8.31 


8.55 


7.94 - 8.52 



a v is the number of binned data points minus the number of fitted parameters. 10 spectroscopic value from Rucinski ct al. (2004), not derived 
photometrically. These values determined using R* = 0.95 ^?Q d This parameter fitted, rather than derived (Rucinski ct al. 2004). 
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TABLE 4 

k 1 Ceti 2004 MCMC fitted parameters 



parameter 


fitted 


p2004 
Parallel 


Parallel Tempering 


p2004 
"■MCMC 


TV /T ^ — i Tl T f 1 






results 


Solution 


D , , 2 
xxeQUCGQ ^ 


n/a 


n/a 


1 Kf) 1 BO 


n/a 


1 46 - 1 52 


i/ a 


n/a 


n/ a 


zoy 


n/a 




1 \ ) 


d&to LlIllCLl 


n/a 


60.6 




60.6 


k 


derived 


n/ a 


n oq^ 111 


n/a 


n 0Q9 n i ni 


vsini (km s — 1 ) b 


derived 


n/ a 


A D(l A oc 

4.80 - 4.86 


n/ a 


4.80 - 4.84 




derived 


n/a 


8.624 - 8.721 


n/a 


8.65 - 8.72 


u 


assumed 


n/ a 


0.6840 


n/a 


0.684 




assumed 


n/a 


0.220 


n/a 


0.220 


<%004 


yes 


1.122 


1.0051 - 1.0150 


1.011 


1.0139 - 1.0150 


E2004.1 (JD-2451545) 


yes 


1.110 


1756.66 - 1756.74 


1.001 


1756.64 - 1756.68 


P2004.1 (days) 


yes 


1.127 


8.704 - 8.773 


1.003 


8.736 - 8.788 


/32004.1 (°) 


yes 


1.183 


9.1 - 22.3 


1.003 


15.6 - 21.7 


72004.1 (°) 


yes 


1.092 


7.51 - 7.80 


1.002 


7.64 - 7.84 


E 20 04_2 (JD-2451545) 


yes 


1.315 


1763.40 - 1763.50 


1.000 


1763.43 - 1763.49 


P2004.2 (days) 


yes 


1.279 


9.050 - 9.162 


1.002 


9.112 - 9.183 


/32004.2 (°) 


yes 


1.162 


-45.8 - -38.6 


1.013 


-47.5 - -43.9 


72004.2 (°) 


yes 


1.074 


11.21 - 15.80 


1.016 


14.65 - 16.83 


E2004.3 (JD-2451545) 


yes 


1.034 


1768.94 - 1769.07 


1.002 


1768.97 - 1769.08 


P2004.3 (days) 


yes 


1.102 


9.518 - 9.598 


1.000 


9.510 - 9.579 


/32004.3 (°) 


yes 


1.203 


71.0 - 77.9 


1.003 


76.6 - 78.0 


72004.3 (°) 


yes 


1.127 


9.04 - 13.11 


1.003 


12.64 - 13.17 


a v is the number of binned data 


, points minus 


the number of fitted parameters. ^ These values determined using i?* 


= 0.95 Rq. 
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Fig. 1. — MOSTlight curves for k 1 Ceti from 2003, 2004, and 2005. The points arc the mean signals from individual 101.413 min satellite 
orbits. The black connected symbols in the lower part of each panel are the the differences between the highest and lowest of 7 adjacent 
simultaneously recorded (background) Fabry images (see text). The 'dip' at 1762 in 2004 corresponds to a total eclipse of the Moon. 
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Fig. 2. — The marginalized likelihood for each of the fitted parameters is shown by the thick curve. The thin dotted curve shows 
the mean likelihood. The best-fitting value (the peak of the distribution) and the 68% credible regions determined from the marginalized 
likelihood are shown by the solid vertical, and dashed vertical lines. The unusual spike in the mean likelihood curve for a low value of k in 
the k histogram is likely due to low numbers of statistics, as indicated by the negligible value of the marginal likelihood curve. 
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Fig. 3. — The solid lines are the 68% and 95% credible regions of the marginalized likelihood for various parameters. The 68% mean 
likelihood credible regions are shown by dotted contours (see text). Obviously there is little correlation between k and i in the bottom 
right panel. 
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Fig. 4. — The best-fitting two spot solution for k 1 Ceti in 2003 (rotating counter-clockwise from top) and seen from the line of sight at 
phases 0.00, 0.25, 0.50, and 0.75 (from left) of the first spot. Middle: the MOST light curve with errors (see text). The continuous line is 
the solution from the "Minimum x 2 " column of Table |31 The dotted line indicates the unspotted normalized signal of the star (U =1.0011). 
Vertical dashed lines indicate phases 0.00, 0.25, 0.50 and 0.75. Bottom: residuals from the model on the same scale. 



14 



Walker et al. 




-0.005 - 
-0.010 - 
-0.015 - 



1750.0 1755.0 1760.0 1765.0 

HJD- 2451 545 



Fig. 5. — The best-fitting three spot solution for re 1 Ceti in 2004 (rotating counter-clockwise from top) and seen from the line of sight at 
phases 0.00, 0.25, 0.50, and 0.75 (from left) of the first spot. Middle: the MOST light curve with errors (see text). The continuous line is 
the solution from the "Minimum x 2 " column of Table |31 The dotted line indicates the unspotted normalized signal of the star (U =1.0149). 
Vertical dashed lines indicate phases 0.00, 0.25, 0.50 and 0.75. Bottom: residuals from the model on the same scale. 
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Fig. 6. — The best-fitting two spot solution for k 1 in Ceti 2005 (rotating counter-clockwise from top) and seen from the line of sight at 
phases 0.00, 0.25, 0.50, and 0.75 (from left) of the first spot. Middle: the MOST light curve with errors (see text). The continuous line is 
the solution from the "Minimum x 2 " column of Table |31 The dotted line indicates the unspotted normalized signal of the star (U =1.0042). 
Vertical dashed lines indicate phases 0.00, 0.25, 0.50 and 0.75. Bottom: residuals from the model on the same scale. 
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Fig. 7. — Top: The 68 % marginalized likelihood contours for the modulus of the latitude, |j8[, and spot period are shown by dots for 
each spot — 2003 (green), 2004 (blue) and 2005 (pink). The central red curve is the solar-period, latitude relation (Equation 1) with the 
mean values taken from Solution 1 of Table 3, k = 0.090 and Peq = 8.77d. The other two red lines correspond to k = 0.085, Peq = 8.74e£, 
and k = 0.096, Peq — 8.81d The vertical dashed lines indicate the range of rotational periods (p = 9.4 ± 0.1 days) found in Ca II H & K 
line reversals over ~30 years. Bottom: The 68 % marginalized likelihood contours for the 2004 data only (cyan) while fitting not for k and 
Peq but explicitly for the periods, P2004_ii P2004_2i and P2004-3- The black curve is the best-fitting curve (k = 0.096 and Peq = 8.68d) 
for the 2004 data, while the red curves are the same as in the upper panel. As can be seen, the black curve fits the data well and thus the 
observed differential rotation pattern is closely similar to solar. 



